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We investigate in detail the possibility of constraining neutrino decays with data from the cosmic 
microwave background radiation (CMBR). Two generic decays are considered vh — > vl4> and vh —* 
vlVlVl- We have solved the momentum dependent Boltzmann equation in order to account for 
possible relativistic decays. Doing this we estimate that any neutrino with mass m > 1 eV decaying 
before the present should be detectable with future CMBR data. Combining this result with other 
results on stable neutrinos, any neutrino mass of the order 1 eV should be detectable. 
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I. INTRODUCTION 

The possibility that one or more neutrino species are 
unstable on a cosmologically interesting timescale has 
been considered very extensively in the literature (see for 
instance PHTl|] and references therein). In particular it 
has been pointed out that cosmology is an excellent lab- 
oratory for neutrino physics and that many exotic neu- 
trino models that are inaccessible to tests in terrestrial 
experiments may be tested using cosmology. 

Big Bang nucleosynthesis has turned out to be a very 
powerful probe for this purpose p-pr|. However, one 
shortcoming of BBN in this regard is that it is only sen- 
sitive to physics before the end of nucleosynthesis, at 
t ~ 10 3 s. After this the relative abundances of light 
elements is fixed until the present. 

Another intriguing possibility is to use the CMBR to 
probe exotic neutrino physics fj^lql. This has been 
done in the past to constrain neutrino decays and other 
effects of non-standard neutrinos. However, there has 
been no real coherent treatment of what limits can be 
put on neutrino decays from CMBR data. In the present 
work we shall try to discuss all the possible decay modes 
of a heavy neutrino in the context of CMBR. 



II. NEUTRINO DECAYS 



vh —* vlvlvl p6| , p7| . Below we shall discuss the possible 
scenarios where such decays can take place. 
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The primary model for this type of decay is the ma- 
joron model |L6j. This is a specific model for the gener- 
ation of neutrino mass. In this model, the neutrino is a 
Majorana particle and <j> is a scalar or pseudoscalar. If 
we assume a Yukawa type interaction of the form 
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and that the daughter particles are massless we obtain a 
rest-frame decay rate of 



r 



9 2 m H 
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It is possible to constrain g from data on neutrinoless 
double beta decays to be g < 10 -4 |fl7|| . Inserting num- 
bers we get a lower bound to the lifetime of 



> 1.6 x 10- ia m^ v s. 



(3) 



This bound is obviously not very restrictive and by far 
the best bounds on such decays come from astrophysical 
considerations. 



A massive neutrino can in principle have many dif- 
ferent decay modes. The simplest possibilities are the 
standard model decays vr - > and vh — > e + e~v e 
[jl(|[l7|. However, they both suffer from the fact that 
their decay products contain electromagnetically inter- 
acting particles and they can therefore be observationally 
constrained. Raffelt [Q provides an excellent review on 
the observational limits on such decays. The main con- 
clusion is that they arc excluded unless the lifetime is 
exceedingly long. 

We are then left with other possibilities which de- 
mand more exotic models. We shall look at two decay 
modes which cover the viable possibilities vh — > v^ab and 



B. VH — > VlVlVl 

This decay is somewhat more complicated. However, 
three-body decays are an equally interesting possibility 
which arises in left-right symmetric models and in gen- 
eral in models with neutral current flavour violation. A 
neutrino decay will then be mediated by either a Higgs, 
Aq, or a Z 1 pq |. In either case, it is safe to assume that 
the boson exchanged is extremely heavy so that the decay 
is effectively a four-point interaction. 

An interesting question is how fast such three-body 
decays can proceed in realistic models. If one assumes a 
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non-derivative coupling as would be the case for both the 
above possibilities then one arrives at an amplitude of 



G 

M = —=U20u\u±Ov3, 
v2 

where the coupling O is given by 



(4) 
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or a combination of these operators. 

Since all of the final-state particles are identical (as- 
suming that there is CP symmetry), all of these differ- 
ent possibilities give the same decay kinematics because 
of the indistinguishability of the final state particles, the 
only difference being the absolute decay rate. This of 
course simplifies the calculations a lot. We need do only 
one calculation for each value of the rest frame decay 
rate, T. This quantity is directly related to the coupling 
strength G by 



i = r = CG 2 ™ 



(6) 



where C is a constant which depends on the specific form 
of the operator O. In the remainder of the paper we shall 
refer to the matrix element for this process as a "weak" 
matrix element because it arises from interactions very 
similar in structure to the standard weak interactions. 

For comparison we also calculate how three-body de- 
cays proceed if the decay matrix element is constant 



E 



Ml 2 = K = G 2 m 4 . 



(7) 



III. FORMALISM 

In the following we will go through the formalism 
needed to calculate decay rates and distributions for both 
two- and three-body decays. In general we need to solve 
the Boltzmann equation which has the form 19 



where 



L[f} = C D [f], 



Tin d f tt df 



(10) 



(ii) 



The right hand side of the equation describes decay terms 
of various sort. These are quite different for two and three 
body decays and we shall discuss each in turn. 



A. two body decays 

Two body decays have quite simple kinematics since 
the rest-frame energy of each daughter particle is well 
defined and equal to mjj/2. We get for the decay terms 



Tm E VH p UH J E 



dE<j,k(f VH , f VL , fy) 
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dE VH A(f UH , f UL , f$) 



This gives a rest frame-decay lifetime of 



(13) 



- = r=^ 2 (8) 

r 51271- 3 

To see what lifetimes can be expected in such a model 
we take an example where the decay is through a flavour 
violating neutral current. We can rewrite the coupling 
constant as G 2 — 9 2 G 2 F and from this we get 

r = 2.88 x 10 14 6> 2 TOm 5 cV s. ( 9 ) 

The limit on is of the order < 10~ 5 |l6| which means 
that for to < 0.2 MeV, the decay lifetime is longer than 
the Hubble time. 

Of course one might postulate a fast invisible decay 
mode to a three neutrino final state, and its effect on 
cosmology would then be one way of detecting it. For this 
reason alone it is interesting to calculate thoroughly how 
such a decay would affect cosmology. In the following 
we shall just take a heuristic approach and examine the 
consequences of a fast three-body decay without looking 
at the theoretical background for this. 



g 777 f "H 

(14) 

where A(f UH , f UL , fy) = f VH (1-/„J(1 + f 4> ) -/^/^(l- 

Ub ) ] m l = m l„- 2 ( m l + m l L ) + ( m l - m l L ) 2 l m ln ■ T is 

the lifetime of the heavy neutrino and g is the statistical 
weight of a given particle. We use g VH = g VL = 2 and 
g$ = 1, corresponding to (f> — <f> The integration limits 
are 

Et H m = ^^[^(l + 4(777 4 /m ) 2 ) 1 /2± (is) 
and 

E t( E » H ) ^ ^[^(l + 4(777 4 /mo) 2 ) 1 / 2 ±p, ff ] (16) 
where the index i — F,<p. 
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B. three body decays 

The decay term in the Boltzmann equation can in this 
case be written as Ji"9 21 



Coif] = ^ / d 3 p 2 d 3 p 3 d 3 p 4 A(f 1 , f 2 , f 3 , U) (17) 
xY,\M\ 2 S\ Pl +p 2 -p 3 - P4 )(2ir)\ 

with d 3 p t = d 3 Pl /(2E i {2n) 3 ) and A(hJ 2 J 3 ,h) = 
(/2/3/41 - /1) - - /a)(l - /a)(l - h))- Details of 
how to evaluate this phase-space integral can for instance 
be found in Rcf. |§. 

The simplest frame to evaluate the phase space integral 
in is the rest frame of the parent particle. Since the 
three body decay scheme does not yield a specific energy 
for the daughter particles, it is interesting to study how 
the energy distribution of daughter particles differ for 
the different interaction matrix elements. In order to 
do this it is very convenient to calculate the differential 
production rate of daughter particles 



d£ 
dp 



= r x 



J 8p/m constant , 

\ 64(p/to) 2 (1 - §p/m) "weak" . 



(18) 



It should be noticed here that the maximum attainable 
momentum for any daughter particle is to/2 in order to 
have 4-momentum conservation. In Fig. 1 we have plot- 
ted the decay distributions for both cases. The most 
interesting feature is that the weak matrix element leads 
to a decay distribution which is stronger peaked in mo- 
mentum space. 
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FIG. 1. The differential energy distribution of the light 
neutrino for three different possibilities. The full line is from 
the "weak" matrix element, the dotted for the constant ma- 
trix element and the dashed for the standard model decay 
h>H — > ULe. + e~ . The lines have been normalised to T = 1. 



In the standard model decay Vh — > v^e^ e~ one gets 
instead the "standard" result for the decay distribution 
of the neutrino 1711 



^- = r x 48(p/to) 2 (1 - \p/m). 
dp 6 



(19) 



The reason why our result is different is that we assume 
that all three final state particles are identical. As for the 
average energy of the emitted neutrinos, it is of course 
equal to m/3 both for the both types of matrix elements. 
For the standard model decay above, it is equal to 7m/20, 
which is slightly higher. 



C. background cosmology 

In order to get the decay into a proper cosmological 
setting it is then necessary to prescribe how the universe 
in general evolves with time. In the simple case where the 
universe is radiation dominated at decay and the energy 
densities of both parent and daughters is subdominant, 
the time-temperature relation is simply pl| 



t(s) = 1.71 g, 



-1/2^-2 

J MeV 



(20) 



where g* is the effective number of degrees of freedom. 
However, in general this is not true because the decay 
might take place while the universe is matter dominated 
or partly matter dominated. In this more general case it 
is necessary to use the equation of energy conservation 



(21) 



which is an exact equation, together with the Friedmann 
equation 



H 2 



inGp k 
~3 W 



(22) 



in order to get a relation between time and temperature 
in the universe Ml. 



IV. NUMERICAL SOLUTION OF THE 
BOLTZMANN EQUATION 

In order to investigate the impact on the CMBR of 
decays it is necessary solve the Boltzmann equation nu- 
merically. To do this we have used a grid in comoving 
momentum space. 

However, first we have to determine the initial condi- 
tions of the system. Primarily this means that we must 
discuss whether or not the daughter particles have ther- 
mal distributions prior to decay. For the ^-particle we 
shall assume that it is not present prior to decay. This 
will be the case if the particle decoupled from thermal 
equilibrium prior to the QCD phase-transition pll. For 
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the light neutrino there are two possibilities, either vl is 
one of the three active species or it is a sterile species. In 
the former case we assume a thermal distribution prior to 
decay and in the latter we assume that there are no light 
neutrinos present before the decay commences. We shall 
always assume that the parent neutrino has a thermal 
population prior to decay. 

For a qualitative discussion of the decays it is conve- 
nient to introduce a relativity parameter a for the de- 
cay Here, we assume that the universe is radiation 
dominated during the decay. A species becomes non- 
relativistic roughly when T = m/3. Using this the rela- 
tivity parameter can be written as 



a = 6.5 x 10 2 gV 2 T s m 



McV- 



(23) 



The decay is then relativistic if a < 1 and non-relativistic 
ifa>l. 

Here we have used the simple relation, Eq. (|2(]). This 
has been done in order to make the decay physics more 
transparent, but in the next section where actual mass- 
lifetime limits are calculated we use the energy conserva- 
tion and Friedmann equations, Eqs. (|2l]-p2|). 
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FIG. 2. The final energy density in decay products for the 
different possible decays, parametrised in equivalent neutrino 
species. The standard model with no decay would correspond 
to JV„, e ft = 1. 



and completely non-relativistic decays. These results are 
compared with the full Boltzmann solution in Fig. 2. The 
full solution clearly has the correct asymptotic behaviour. 



B. three body decays 



For non-relativistic decays, the cosmic frame is equal 
to the rest-frame of the parent particle. In this case, for 
equal total decay rates, the "weak" and the constant ma- 
trix element give somewhat different decay distributions. 
The reason is that energy and number conservation only 
yields two fixed parameters. For an equilibrium distribu- 
tion this is sufficient to describe the distribution fully by 
a temperature and a pseudo-chemical potential, but the 
distribution of the daughter particle is not an equilibrium 
one so that there are no two parameters fully describing 
the distribution function. 

For extremely relativistic decays the decay and inverse 
decay installs equilibrium in the particle distributions for 
both parent and daughter. In this case, the decay pro- 
ceeds in complete equilibrium and the final distribution 
functions become identical, independent of the actual 
matrix element. 

We have plotted the final momentum distribution of 
the daughter particles for different values of the relativ- 
ity parameter, a. Indeed one sees that for very relativis- 
tic decays the distributions become identical. For non- 
relativistic decays they become very significantly differ- 
ent, even though the zeroth and first moments are iden- 
tical. If one is interested in applications where the actual 
shape of the distribution is important, as would be the 
case if the daughter is an electron neutrino present dur- 
ing nucleosynthesis [|oJTJ] or an eV particle constituting 
the Hot Dark Matter [ 23| . this difference in distribution 
could potentially be important. However, if one is inter- 
ested only in the final number and energy density, the 
difference between using different matrix elements is in- 
deed very small and can be safely ignored. 

It is quite interesting to compare the final decay dis- 
tributions with equilibrium distributions of the form 



/eqO) 



1 



ex P 



(24) 



A. two body decays 

The two-body decay is interesting in that one of the 
final state particles is a boson. This means that because 
of stimulated emission effects there will be a large popu- 
lation of low momentum 0's. This of course only happens 
if these states are accessible, i.e. for relativistic decays 
f|po"[ In the Appendix we have calculated the final en- 
ergy density in decay products for both very relativistic 



with the same number and energy densities as the true 
distributions. We have plotted these equilibrium dis- 
tributions together with the full distributions calculated 
from the Boltzmann equation in Fig. 3. It is seen that 
in the strongly relativistic case the decay proceeds in 
equilibrium and the final distributions are of equilibrium 
form. For strongly non-relativistic decays the equilib- 
rium distributions are quite poor approximations to the 
true distribution functions, as expected. 
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FIG. 3. The final momentum distribution of the light neu- 
trino after decay for the three body decay vh — » vlVlVl- The 
full line is for the "weak" matrix element, the dotted for the 
constant matrix element and the dashed line shows an equilib- 
rium distribution with the same number and energy density. 
The five sets of curves correspond from left to right to a — 
0.2, 2, 10, 100, 500. 

V. CMBR EFFECTS 

Using the results on the energy density evolution of the 
decaying particle and its daughters we can now proceed 
to investigate the CMBR effects from neutrino decays. 
The fluctuations are usually described in terms of spher- 
ical harmonics [03 

T(9,(f>)=J2 a lmYi m (0,4>), (25) 
lm 

where the coefficients are related to C\ coefficients by 

Q = (\a lm \ 2 ). (26) 

These fluctuations were first detected in 1992 by the 
COBE satellite ||, but only for I < 20. At such low 
I the power spectrum is almost degenerate in the cos- 
mological parameters and no real constraints are obtain- 
able. In the next few years, however, the power spectrum 
will hopefully be measured out to I ~ 2500 by two new 
probes, MAP and PLANCK |||, and using this data 
should yield precision measurements of the physical pa- 
rameters at recombination. We shall here assume that 
MAP will deliver data out to Z = 1000 and PLANCK to 
I = 2500. 

Since we do not yet have the precision data in hand it 
is hard to foretell how accurately the different parame- 
ters can actually be measured. The usual procedure is to 
use what is called error forecasting |27j . This method as- 
sumes an underlying cosmological model to be the "true" 



model. From that one can calculate how sensitive the 
data are to changes in the cosmological parameters. It 
should be noted here that this method can at best give 
an estimate of the obtainable precision. Once the data 
becomes available the actual problem of determining the 
cosmological parameters will be much tougher because 
the whole parameter space has to be investigated. 

In the present paper we shall assume as the reference 
the standard cold dark matter model which can be de- 
scribed by a vector of values in the many-dimensional 
parameter space |27j 

e = (fi,fi fe ,A,/i,n,r). (27) 

As free parameters we have chosen the total density, fi, 
the density in baryons, fib, the cosmological constant, 
A, the Hubble parameter, h, the spectral index, n, and 
the optical depth to reionisation, t. In addition to this 
we have included as free parameters the energy density 
injected by neutrino decays as well as the time at which 
the energy is injected. 

Our choice of reference model is "standard" cold dark 
matter 

Ocdm = (1,0.08,0,0.5,1,0). (28) 

Note that we have chosen a limited set of free parameters. 
In principle one should use a parameter space consisting 
of all possible free parameters. 

Our procedure will then be the following: for given life- 
time th of the heavy neutrino we calculate which mass 
range will be detectable by the upcoming CMBR exper- 
iments. To estimate the obtainable precision we use the 
Fisher information matrix [£7j 

la = f> + 1) [Ci + a. mm ]- 2 (29) 

where C/ i0 rror represents the experimental error which we 
neglect in the present paper. The standard error in pa- 
rameter i is then 

°i * (1-%, (30) 

if all cosmological parameters must be determined simul- 
taneously, and 

° 2 i-{hi)-\ (31) 

if all the other parameters are assumed to have already 
been determined |l2| . To calculate actual CMBR spectra 
we use the CMBFAST package designed by Seljak and 
Zaldarriaga pq ). 

To be on the conservative side we have calculated the 
sensitivity of the CMBR for the case where all parame- 
ters must be determined simultaneously. The results are 
shown in Fig. 4 where the upper panel shows the expected 
accuracy of PLANCK and the lower panel that of MAP. 
For very small neutrino masses the the decay cannot take 
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place before t(m/3). Also, the CMBR data are less sensi- 
tive to late decays. Therefore there is a sharp and well de- 
fined minimum mass which is detectable through CMBR 
measurements. This lower limit depends on the decay 
mode because the energy deposited in relativistic parti- 
cles depends on the specific decay kinematics. 



data, which is of low accuracy. What we have shown 
is that the CMBR is a very sensitive probe of neutrino 
decays. Clearly, our limits only apply to the case where 
decay takes place after neutrinos have decoupled com- 
pletely from thermal equilibrium at T ~ 0.1 MeV. 
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FIG. 4. The expected detectable neutrino mass and life- 
time with new CMBR data. The region to the right of the 
curves is the detectable range. Curve labels are as in Fig. 2. 
The shaded regions in both plots correspond to the detectable 
range found by Hu, Eisenstein and Tegmark pS[ . 

For both early and late decays the main effect on 
CMBR is through the ISW effect @- This effect is due 
to the fact that if the universe has a significant radiation 
content the linear gravitational potential is not constant 
in time. This leads to an enhancement of the CMBR 
fluctuations, primarily on the scale corresponding to the 
horison size when the energy is injected. 

Our results may be combined with the results of Hu, 
Eisenstein and Tegmark |2g| on stable low mass neutri- 
nos. They estimate that for stable neutrinos a combina- 
tion of large scale structure surveys and CMBR measure- 
ments should be able to measure a neutrino mass of 1-2 
eV. 

Our estimate of the obtainable precision is clearly 
much more restrictive than that found by Lopez et al. 
[ ]l"3| ], because they rely only on already existing CMBR 



VI. DISCUSSION 

In this paper we have calculated the constraints that it 
should be possible to put on neutrino decays using data 
from the new CMBR experiments. We have paid par- 
ticular attention to the different possible decay modes, 
discussing both two- and three-body decays in detail. 
Also, the full Boltzmann equation was solved, instead 
of the momentum integrated version [5jJ . This has the 
clear advantage that it is possible to do calculations for 
relativistic neutrino decays where inverse decays are in- 
cluded. 

Using this formalism we estimated that it should be 
possible to constrain any type of neutrino decay if the 
mass is of the order 1-2 eV. If the lifetime is short it may 
be possible to detect neutrino masses down to 0.1 eV 
using this method. Our results complement and extend 
existing results on CMBR constraints on neutrino decays 

Finally, it is of interest to discuss how a decaying neu- 
trino scenario relates to the recent results from Super- 
Kamiokande. According to these results the muon neu- 
trino oscillates with near maximum mixing with another 
neutrino, most likely the tau neutrino [^0|. The mass 
difference between these two states is of the order 0.1 
eV, and if one simultaneously believes the solar neutrino 
anomaly to be due to v e — oscillations |3l| there is no 
room left for cosmologically interesting neutrino decays. 
However, there might be other possibilities. For instance 
the solar neutrino solution might be due to oscillations 
with a sterile neutrino state. The only constraint is then 
that f M and v T should be almost degenerate in mass. 
This could still allow for cosmologically interesting and 
detectable decays. Until there are more firm results on 
this subject, neutrino decays remain a viable possibility. 



APPENDIX A: ASYMPTOTIC BEHAVIOUR 

1. decay in equilibrium 

If the decay is relativistic it is a very good approxi- 
mation to assume that it proceeds in equilibrium. This 
means that the distribution functions are all of equilib- 
rium form 



feq,i(E) — 
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In this case, instead of using the Boltzmann equation, we 
can derive a very simple set of equations for the temper- 
ature and pseudo-chemical potentials. 

The procedure is the following: One starts at some 
initial time with given initial conditions. Then at each 
timestep the distributions are evolved forwards in time 
and then allowed to equilibrate completely before the 
next timestep is taken. 

We start with the two-body decay. In this case the 
equations to be solved at each timestep are 



m = Pl (A2) 

no = n H + n L (A3) 

Pa = Ph + PL + P<t> (A4) 

n L = n<f,. (A5) 

T H = T L = (A6) 



An interesting possibility in this case is that it is possible 
to form a Bose-condensate from the pseudo-scalar parti- 
cle. This phenomenon has been investigated previously 
in the literature || 

For the three body decay, the corresponding equations 
are 



2. non-relativistic decays 

For strongly non-relativistic decays the behaviour is 
also quite simple. Since the energy in the decaying par- 
ticle is equal to the mass, m, we can write an equation 
for the evolution of energy density in decay products 

-=f— (AH) 
Pa J ti T ( E o) 

where po is the energy density in a standard massless 
species. (Eo) = 3.151 T is the mean energy of a massless 
fermion and ti is some initial time which is taken to be 
ti <C t. Solving this equation we get 

_p = 270C(3) al/2 ^ _ 84al/2 (A12) 

PO llt'l 1 

for decays in the radiation dominated era and 

^ = 540C|3) r /5V 2 /3^ . 86a2/3; (A13) 
Po 7 7r \3J 

for decays in the matter dominated regime. 



p H = iph 
n = n H + n L /3 
Po = Ph + PL 

T H = T L . 



(A7) 
(A8) 
(A9) 
(A10) 



We have solved these equations to obtain the final en- 
ergy density in the daughter products since this is the 
relevant quantity to know if we are interested in effects 
on the CMBR. 

It is of course also worth noting that for relativistic 
decays there is no entropy production since all the distri- 
bution functions stay in kinetic equilibrium throughout 
p9[ . Also, the final energy density in decay products 
is completely independent of the expansion rate of the 
universe since equilibration is assumed to happen instan- 
taneously at each temperature. In Table I we have shown 
the final energy density for relativistic decays for the four 
different possibilities which we have treated. 

TABLE I. Final energy density in the decay products in 
units of a massless neutrino species. 



Decay 






vh —* vl4> 





1.20 


VH — > VL<t> 


thermal 


1.48 


VH — » VlVlVl 





2.20 


vh — » vlvlvl 


thermal 


2.63 
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